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We combine the BCS self-consistency condition, a semiclassical expansion for the spectral density 
and interaction matrix elements to describe analytically how the superconducting gap depends on 
the size and shape of a 2d and 3d superconducting grain. In chaotic grains mesoscopic fluctuations 
of the matrix elements lead to a smooth dependence of the order parameter on the excitation energy. 
In the integrable case we find shell effects i. e. for certain values of the electron number N a small 
change in N leads to large changes in the energy gap. With regard to possible experimental tests we 
provide a detailed analysis of the dependence of the gap on the coherence length and the robustness 
of shell effects under small geometrical deformations. 
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Finite size effects are well documented [l| in fermionic 
interacting systems such as atomic nuclei and atomic 
clusters. It is also well established 0, Q that the more 
symmetric the system is, the stronger are these correc- 
tions. For instance, the existence of magic numbers sig- 
naling the presence of a particularly stable nucleus has its 
origin in the gap between the ground state and the first 
excited states caused by the high degree of symmetry of 
the system. 

In the field of mesoscopic superconductivity, the study 
of finite size effects also has a long history. Already 
fifty years ago, Anderson noted [3| that superconductiv- 
ity should break down in small metallic grains when the 
single particle level spacing at the Fermi energy is com- 
parable to the bulk superconducting gap. In the sixties 
the size dependence of the critical temperature and the 
superconducting gap were studied in for a rectangular 
grain in [5[ and for a nanoslab in|6j. Thermodynamical 
properties of superconducting grains were investigated 
in [3]. Results of these papers are restricted to rectan- 
gular grains, and superconductivity is described by the 
Bardeen, Cooper, and Schriffer (BCS) theory 

The experiments by Ralph, Black, and Tinkham in 
the mid nineties [9( on Al nanograins of typical size 
L ~ 3 — 13 nm showed that the excitation gap is sen- 
sitive to even-odd effects. More recently it has been ob- 
served [l(| that the critical temperature of superconduct- 
ing ultra-thin lead films oscillates when the film thickness 
is slightly increased. These results have further stim- 



ulated the interest in ultrasmall superconductors (Tll - 
[bfij . For instance, pairing, not necessarily BCS, in a har- 
monic oscillator potential was investigated in 13]. The 
critical temperature and the superconducting gap for a 
nanowire were reported in [l4j by solving numerically 
the Bogoliubov - de Gennes equations. In [15( the super- 
conducting gap and low energy excitation energies in a 
rectangular grain were computed numerically within the 
Richardson model [lq |. Shell effects in superconducting 
grains with radial symmetry were studied theoretically 
in [13, EH HH ■ Recent experimental results [l9[ in semi- 
spherical Sn nanograins have confirmed that shell effects 
induce strong deviations in the energy gap with respect 
to the bulk limit. Strong fluctuations of the energy gap 
as a function of the system size have been observed with 
a maximum enhancement of about 60% for sizes ~ lOnm. 
Mesoscopic corrections to the BCS energy gap were also 
considered in (2^, [2l| ■ 

We note that if the mean single particle level spacing is 
larger than the bulk superconducting gap, the BCS for- 
malism breaks down. However, an analytical treatment 
is still possible [12] with the helpof an exactly solvable 
model introduced by Richardson JJ3J in the context of nu- 
clear physics. In particular, finite-size corrections to the 
predictions of the BCS theory have been recently studied 
in US-El. 



Despite this progress, a theory that accounts for all 
relevant mesoscopic effects in superconducting grains has 
not emerged so far. The Richardson model alone cannot 
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provide the foundation for such a theory as it does not 
allow for mesoscopic spatial fluctuations of the single par- 
ticle states. In the present paper, for the particular cases 
of chaotic and rectangular shaped grains, we develop such 
a theory based on the BCS theory and semiclassical tech- 
niques. This formalism permits a systematic analytical 
evaluation of the low energy spectral properties of super- 
conducting nanograins in terms of their size and shape. 
Leading finite size corrections to the BCS mean field can 
also be taken into account in our approach, see [28| for 
further details. Results for 3d grains were also previously 
published in Q . Here we discuss both the 2d and 3d 
cases as well as provide a more detailed account of the 
techniques utilized. Moreover, we study the dependence 
of the mesoscopic BCS order parameter (superconduct- 
ing gap) on the coherence length, and the robustness of 
shell effects. 

For chaotic grains, we show that the order parameter 
is a universal function of the single particle energy, i.e. 
it is independent of the particular details of the grain. 
The mesoscopic fluctuations of the matrix elements of 
the two-body interactions between single particle eigen- 
states are responsible for most of the deviations from the 
bulk limit. For integrable grains, we find that the su- 
perconducting gap is strongly sensitive to shell effects. 
Namely, a small modification of the grain size or num- 
ber of electrons inside can substantially affect its value. 
Throughout the paper we study clean (ballistic) grains. 
The mean field potential is thus an infinite well of the 
form of the grain. We restrict ourselves to system sizes 
such that the mean level spacing around the Fermi energy 
is smaller than the bulk gap, so that the BCS formalism 
is still a good approximation. For the superconducting 
Al grains studied by Tinkham and coworkers @, this 
corresponds to sizes L > 5 nm. 

Our results are therefore valid in the region, kpL ^> 1 
(limit of validity of the semiclassical approximation [U, 
Q), 5/A < 1 (limit of validity of the BCS theory), and 
( > ( > I> (condition of quantum coherence). Here kp, 
£ = Hvp/Aq, /, S, An are the Fermi wave vector, the 
superconducting coherence length, the coherence length 
of the single particle problem, the average single particle 
level spacing, and the bulk gap. The Fermi velocity is 
vp = hkp/m. Conditions kpL 3> 1 and 5 / Aq < 1 hold 
for Al grains of size L > 5 nm. Further, in Al grains 
£_~ 1600 nm and I > 10 4 nm at temperatures T < 4K 
24]. Therefore, the above region is well accessible to 
experiments. 



I. THE SUPERCONDUCTING GAP IN THE 
BCS THEORY 



where c na annihilates an electron of spin a in state n, 

I n>n , = I(e n ,e n ,) = XV 6 J ^ 2 n (r)ip 2 n ,(f)dr (1) 

are matrix elements of a short-range electron-electron in- 
teraction, A is the BCS coupling constant. ip n and e„ are 
the eigenstates and eigenvalues of a free particle of ef- 
fective mass to in a clean grain of volume (area) V (A). 
Eigenvalues e n are measured from the actual Fermi en- 
ergy tp of the system. In this notation the mean level 
spacing is S = 1/^tf(0), where ^tf(O) is the spectral 
density at the Fermi energy in the Thomas-Fermi ap- 
proximation. 

The BCS order parameter is defined as 



A(e r , 



Within BCS theory, it is determined by the following 
self-consistency equation [25j: 



E 

n >\<e D 



(2) 



where en is the Debye energy. This result is obtained in 
the grand canonical approximation [8j. Note that, the 
BCS order parameter A n is an explicit function of the 
single-particle energy e n since the matrix elements /(e, e') 
are energy dependent. 

Introducing the exact density of single-particle states 
2/(e') = J2n' <K e ' — e «')i one can wr ite Eq. ([2]) in integral 
form, 



A(e) 



A(e')/(e,e') 



A 2 (e') 



v{e )de . 



(3) 



The gap equation ([3]) will be the main subject of our 
interest. As soon as the order parameter A(e) is known, 
the low lying (single-particle) excitation spectrum, E = 
\J A(e) 2 + e 2 , is also determined. 

In the large volume (area) limit, the spectral density, to 
leading order, is given by the Thomas-Fermi expression 



^TF(e') = 2 




for 3d 
for 2d, 



(4) 



where the factor two in front stands for spin degeneracy. 
In addition, in the bulk limit the matrix elements ([1]) for 
chaotic grains are simply I(e, e') = XS as a consequence of 
quantum ergodicity. The gap is then energy independent 
A(e) = A , and Eq. © yields the BCS bulk result, 



An 



2e D e~ 



(5) 



Throughout the paper pairing between electrons is de- 
scribed by the BCS Hamiltonian, 



H — ^ e n c' na c rl 



^i,n' c it c n4. c "'4- c n't' 



As the volume of the grain decreases, both v(e') and 
/(e, e') deviate from the bulk limit. In this region a more 
general approach to solve Eq. ((3]) is needed. 

Since we are interested in the regime of many particles 
(z/tf(0)cf S> 1), an appropriate tool is the semiclassical 
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approximation in general and periodic orbit theory [29j in 
particular (see the Appendix for an introduction) . These 
techniques yield closed expressions for v(e') and J(e, e') 
in terms of quantities from the classical dynamics of the 
system, which allows us to calculate analytically the re- 
sulting superconducting gap. Such explicit expressions 
for the superconducting gap enable us to study devia- 
tions from the BCS theory, the spatial dependence of the 
gap, and the relevance of shell effects in realistic, not 
perfectly symmetric grains. 

Our general strategy can be summarized as follows: 

1. Use scmiclassical techniques to compute the spec- 
tral density v(e') — J2 n ' ^( e ' " e ™') anc ^ -^( e > £ ') as 
series in the small parameter 1/kpL, where Uf is 
the Fermi wave- vector and L ~ l/ 1 / 3 (~ A 1 / 2 ) is the 
linear size of the grain (section [TTI and Appendix). 

2. Solve the BCS gap equation ([2|) order by order in 
1/kpL (Section |in|. 

3. Study the impact of small deformations of the 
shape of a symmetric grain on the gap in realis- 
tic models of the grain (Section IIV[) . 

Finally we stress that all the parameters in our model 
X,kp,eD,£F are the actual parameters that characterize 
the material at a given grain size and not necessarily the 
ones at the bulk limit. 



for Dirichlet (— ) or Neumann (+) boundary conditions. 
In Eq. ([TJ, S is the surface area of the 3d cavity and C 
its mean curvature, while C is the perimeter in the 2d 
case. 

The oscillatory contribution to the density of states is 
given by the Gutzwiller trace formula [29$ . 

i-^-y"' A e l ^ kFL f + ^e i ^ kFLp 3d 
T « , , (8) 
T^E p M l[kFLp+M ^ F p 2d. 

The summation over classical periodic orbits (p) with 
length L p only includes orbits shorter than the quantum 
coherence length I of the single-particle problem. The 
scmiclassical amplitude A p and phase j3 p in Eq. ([5]) can 
also be computed explicitly using the knowledge of peri- 
odic orbits. As was mentioned previously the parameters 
Uf and cf in the above expressions refer to the Fermi 
wavevector and Fermi energy of the system at a given 
grain size. Within the free Fermi gas approximation it is 
possible to relate the bulk Fermi energy with the one at 
a given finite size by simply inverting the relation 



v(e)de 



(9) 



where v(e) is the spectral density and N is the number 
of particles. 



II. SEMICLASSICAL APPROXIMATION FOR 
THE DENSITY OF STATES AND INTERACTION 
MATRIX ELEMENTS. 

The first step to solve the gap equation is to find ex- 
plicit expressions for the spectral density v(e') and the 
interaction matrix elements I(e, e') as series in a small 
parameter 1/kpL. While the semiclassical approxima- 
tion for the spectral density has been known for a long 
time [29(, the calculation for the matrix elements has 
only recently attracted some attention [28], [3(| ■ Here we 
state the results and refer the reader to the Appendix for 
details. 



B. Matrix elements 

The calculation of the interaction matrix elements 
/(e, e') is more complicated as it requires information 
about classical dynamics beyond periodic orbits. For a 
chaotic cavity the final result (Appendix l"2l . 



■ 1 + 



(10) 



rshort 



(0) 



16*4 v 2 



l + I|h°*(0,e-eO+4? S (0^-e') 



3d, 



2d, 



has two types of contributions. Identical pairs of short 
classical trajectories hitting the boundary once give 



A. Spectral density 

In the semiclassical approximation (see Appendix [TJ , 
the spectral density is given by 



Ke')^TF(O) [l + g(0) + W)]> 



(6) 



with a monotonous g(e') and oscillatory tjf(e') (as func- 
tions of system size) parts. The notation g(e = 0) means 
that g is evaluated at the Fermi energy. This contribution 
is given by the Weyl expansion [l(, 



9(0) = 



±— 

. 2k F A' 



3d, 
2d, 



(7) 
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ik F V 



^ ort (0, £ -e') = ^ 



27rk F A 



c 



i(lfc|^)-Ci( 



qi . Si(4k F L 



3d, 
2d, 



(11) 



with C = 0.339... a numerical constant given in the Ap- 
pendix, and Ci(x) the cosine-integral function. 

In the so-called diagonal approximation (see Ap- 
pendixl"2l the contribution of longer classical trajectories 
is 



flong/ 
7 d e \ € F,£ 



in,(^) 3d, 




(12) 
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where 

Ui(w) = f^2D^cos[wk F L^(r)]dr (13) 
■* 7 

is an integrated sum over trajectories 7(r) starting and 
ending at position r. As detailed in Appendix \2\ due 
to the ergodicity of the chaotic classical systems, in the 
limit I S> L, Eq. (fT3"j) simplifies to 

sin (wkpl) 

ig; — s — 

For integrable grains there is no universal expression 
for /(e, e'). We restrict ourselves to the rectangular geom- 
etry where to a good approximation the matrix elements 
are energy independent. 

Using the knowledge of v(e') and I(e,e') as series in 
1/kpL, we solve the gap equation @ in different situa- 
tions of interest. The resulting gap function, in general, 
depends the single-particle energy e, the size of the sys- 
tem, and the number of particles (or, equivalently, Fermi 
energy e F ). 

III. SOLUTION OF THE GAP EQUATION IN 
THE SEMICLASSICAL REGIME 

In this section we solve the gap equation Eq. ([3]) for 
A(e). For a rectangular box in two and three dimensions 
the gap equation is algebraic, since A(e) = A is energy 
independent. In the chaotic case, however, we get an 
integral equation due to the energy dependence of the 
interaction matrix elements. As we will see, both cases 
can be solved analytically order by order in 1/kpL. 

A. Rectangular box in two and three dimensions 

For the rectangular box the matrix elements are 

J(e,e')= II {l + 5 HA /2)/V (15) 

i=x,y,z 

where ej oc kf, pi — hki is the conserved momentum in 
the i — x,y, z direction and here S stands for Kronecker's 
function. We first investigate the role of these matrix ele- 
ments on the energy gap. Qualitatively we expect an en- 
hancement as I{e,e') > 1/V. This enhancement should 
not be large for S/Aq < 1 as the spectrum of a rectangu- 
lar grain has only accidental degeneracy, namely, ti — d i 
typically implies that i — i' . For a perfectly cubic grain 
the enhancement is expected to be larger due to level 
degeneracy although they will still relatively small for 
<5/A <C 1. The numerical results of Fig. Q] (upper) for 
the gap as a function of the grain size confirm this pre- 
diction. We compare the cases of trivial matrix elements 
I(e, e') « 1/V and Eq.(fT5]) (see caption for details). In 



2nr 




L(nm) 

FIG. 1. Upper figure: The energy gap A in units of the bulk 
gap A for a cubic grain of side L with A = 0.3, eo = 32meV, 
£f ~ 11.65eV, &f = 17.5nm _1 as a function of the grain size 
L. The chemical potential was computed exactly as a function 
of TV by inverting the relation — v(e)de where u(e) is 
the spectral density. Similar results (not shown) are obtained 
for other values of A. Red circles stand for the exact numeri- 
cal solution of the gap equation Eq.([2} with matrix elements 
Eq, (|15[l . The black curve is its average value A ave . Blue 
squares is the numerical solution of Eq.Q for trivial matrix 
elements I(e,e') — 1/V. The green curve is its average value. 
Lower figure: we represent the standard deviation of the gap 
cr(L) in units of the average gap A at)e , a typical estimation of 
the average fluctuation, as a function of the grain size. The 
black (green) curve is the typical deviation for the case of non 
trivial matrix elements given by Eq.{TS]) (J(e,e') = 1/V). As 
can be observed, in the region S/Aq S> 1 (Al L S> 6nm), 
in which our semiclassical formalism is applicable, the non- 
trivial matrix element Eq. (|15p does not modify substantially 
the average gap or the typical fluctuation. We note that the 
average fluctuation (see also figure 2) is in reasonable agree- 
ment with the theoretical prediction, m J j^- J2jj. 



the region in which our results are applicable 6 <C Ao 
(L ^S> 6nm) the enhancement of both the gap average 
(upper plot) and fluctuations (lower plot) due to Eq. (fT5| 
is small. We note that in the numerical calculation the 
chemical potential is not the bulk Fermi energy but it is 
computed exactly for each grain size (see caption). This 
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induces an additional enhancement of the average gap 
with respect to the bulk limit A . 

Since we are mainly interested in the study of gap fluc- 
tuations (see below) we neglect in the rest of this section 
the non trivial part of Eq.(fT5]) (J(e,e') « l/V). There- 
fore to a good approximation the gap does not depend 
on energy, A(e) = A, and satisfies the equation, 



Ve' 2 + A 2 



(16) 



where g(0) for a 3d rectangular box is given by Eq. ([7]) 
without the curvature term. 

Using Eq. (0 for g(0) and Eq. (JS) for gi(e') (from 
now on we drop the subscript I to simplify the notation), 
and taking into account the scaling of each contribution 
with 1/fcpL as described in the Appendix, we look for a 
solution of the gap equation (fl6|) for the 3d case in the 
following form: 



A = A (l + /« + /^ + /^), 



(3/2) 



f(2h 



(17) 



where /W oc \/{k F L) n . Substituting A into Eq. CT . 
expanding in powers of l/Zc^L, and equating the coeffi- 
cients at each power, we obtain an explicit expression for 



A/ 



(i) 



7(0) + 



(3/2) _ 



E~ 



--D o (2) (e') 

de', 



-.de' 



s D Ve' 2 + Aq 



(18) 



A/ 



(2) _ 



? 2 y_ £D v^Ta 2 

(/ (1) -.9(0) 

V 2 1 £D ( £ ' 2 + A 2 )3/2 



(19) 



l2 reD -(i) (c/) 



de', 



where oc (kpL)~ k denotes the oscillating part of the 
spectral density. Explicit expressions for g^ k \ g\ , and 
gfj for a rectangular box in terms of periodic orbits can 
be found in the Appendix and also in Ref. [l|. 

Equations (fT8|) and (fT9|) can be further simplified by 
the following argument. After we express g^ 3 ',g^ 2 ' and 
§W in terms of a sum over periodic orbits, the integration 
over e' can be explicitly performed. The resulting expres- 
sion is again an expansion in terms of periodic orbits with 
two peculiarities: a) the spectral density is evaluated at 
the Fermi energy and b) in the limit ejj 3> Ao the con- 
tribution of an orbit of period L p is weighted with the 
function 



A f c 
W(L p /0 = -J 



cos{L p t/£) 



dt. 



(20) 



This cutoff function is characteristic of the BCS theory as 
opposed to the smoothing due to temperature or inelastic 
scattering (recall that in this paper we assume that the 
single-particle coherence length I is much larger than su- 
perconducting coherence length £). In a similar fashion, 
the last term in is weighted with 



W 3/2 (L p /0 



A 2 



cos(L p t/£) 

-co (1+i 2 ) 3 / 2 



dt. 



The effect of W 3 /2{L p /^) is, again, to exponentially sup- 
press the contribution of periodic orbits longer than £. 
Therefore the sum over periodic orbits in the definition 
of the spectral density is effectively restricted to orbits 
with lengths of the order or smaller than the supercon- 
ducting coherence length £. 

Following standard semiclassical approximations, we 
introduce 5^(0) as a spectral density evaluated at the 
Fermi energy with a cutoff function that suppresses the 
contribution of orbits of length L p > £ . With these defi- 
nitions, we get 



A/ (1) = [s(0)+4 3) (0) 

^ (3/2) = E.9X(o), 

A/^E^ ) 



(21) 



+ / 



(i) 



/W- g(0) -E^O) 



Eq. ([21]) is our final result for the finite size corrections to 
the gap function for a 3d rectangular box. As expected, 
it is expressed in terms of classical quantities such as the 
volume, surface, and periodic orbits of the grain. 

In Fig. [5] we compare the analytical expression for 
the gap (TrT)) and (|2ip (solid blue line) to the numerical 
solution of the gap equation using the exact one-body 
spectrum (circles) and the semiclassical prediction for 
the spectral density (red squares). It is observed that 
the analytical expression for the gap is in fair agreement 
with the exact numerical results. Moreover it is also clear 
from the figure that the semiclassical formalism provides 
an excellent description of the numerical results. We note 
that the small differences observed for small values of the 
gap are a consequence of the finite I ~ 50R single parti- 
cle coherence length entering in the semiclassical expres- 
sion of the spectral density Eq.([5]). Since our motivation 
here is test the validity of the semiclassical formalism we 
are assuming for simplicity that the chemical potential is 
fixed at the bulk Fermi energy. 

The following argument can shed light on our results. 
The density of states cannot be pulled out of the energy 
integration in the gap equation (|16p unless it is smoothed. 
However, this is exactly what our result Eq. (1211) means, 
since truncating the sums is equivalent to smoothing the 
energy dependence. We conclude that our result Eq. (CR 
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FIG. 2. The energy gap A in units of the bulk gap Ao for 
a cubic grain with A = 0.3, en = 32meV, ~ 11.65eV, 
kp = 17.5nm -1 as a function of the number of particles N 
(L ~ 13.23 — 13.32nm) inside of the grain. The solid line is 
the analytical prediction from (|17|l and (I21|l . Black circles ( 
red squares) are results from a numerical evaluation of the 
gap equation using the exact (semiclassical Eq. (J5|) spectral 
density. The semiclassical formalism provides an excellent 
description of the exact numerical results for the gap. We 
stress that, for the sake of simplicity, it has been assumed 
that/ = l/V. 



should be similar to the standard BCS solution in the 
bulk, A = 2e iD e~ 1 / A , with the substitution A — > A(l + 
g(0) + <?£(0)). Indeed, an expansion of this expression in 
1/kpL gives exactly Eq. (f2"TT) . 



In order to simplify notation from now on we will drop 
the subscript £ in the spectral density <?{ smoothed by 
the cutoff function W(L p /£). In 2d we find, 



A = A (l + /( 1 / 2 > + /W), 



(22) 



with 



A/^)=^)( ) 
A/ (1) -5(0)+^3, (1) (0) 



(23) 



z=l,2 

1 A • • ■ i - 



A 



[ffS(0) 



The sums implicit in gi,(ji,j are smoothly truncated by 
the same weight function W(L p /£). Similar to the 3d 
case, the above result can also be obtained by expanding 
the bulk expression for the gap with the full density of 
states in l/(kpL). We note that, contrary to the 3d case, 
in 2d grains, oscillatory contributions to the density of 
states are of leading order. 



B. 3d chaotic cavity 

The energy dependence of the interaction matrix ele- 
ments, J(e, e'), in this case is given by Eqs. (p!T71 FT4|) . i.e. 



where 



ttS 



7T 2 S 2 



4k F V 16k 2 F V 2 V 
4-7T 2 sm(kFloj) 



1 



L.3 
ft p 



(24) 



The details of the calculation based on the semiclassi- 
cal approximation for Green's functions can be found in 
Appendix [2J 

The above expression for I(e, e') together with the 
semiclassical expression for the spectral density © are 
the starting point for the calculation of the supercon- 
ducting order parameter. The energy dependence of the 
matrix elements implies a gap equation of integral type 
and, most importantly, that the order parameter itself 
depends on the energy. Based on the 1/Hf L dependence 
of the different contributions to I(e, e'), we write 



A(e) = Ao l + / (1) +/ w +/ w (e) 



f(2) 



f(3), 



(25) 



for a 3d chaotic grain. Substituting this expression into 
the gap equation ((3|) and comparing powers of 1/kpL, we 
get a simple algebraic equation for f^ 1 ' with the solution 



(i) _ 



(1±1) 



Stt 
4k f V 



(26) 



It shows that for Dirichlet (-) boundary conditions, the 
superconducting order parameter for a chaotic 3d cav- 
ity does not have mesoscopic deviations of order 1/fcp-L. 
This suppression is a hallmark of the chaotic case and 
appears due to the fluctuations of the interaction ma- 
trix elements. It can be also found by substituting 
A — » A(l + STr/AkpV) into Eq. ([5]), which accounts only 
for the surface contribution to the density of states, and 
expanding the modified Ao to first order in 1/kpL [21] . 
The second order correction reads 



with 



.9(0) 



k 2 F V 



2n 



A J \Ak F V 



+ 5(0), (27) 



J2A p W(L p /$cos(k F L p +p p ), (28) 



where the contribution of periodic orbits L p longer than 
the coherence length £ is exponentially suppressed. 

Equating terms of order (kpL)~ 3 , we obtain for 
f^ 3 \e) an integral equation of the form /®(e) = 
h(e) + J K(e')f^(e')de' , which is solved with the ansatz 
j( 3 )(e) = h(e) + c, where c is a constant. We obtain 



f i3) (e) 



ttXS 
Ao" 



Ao 



A 2 



(29) 
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Note that a) since 5/Aq < 1 is an additional small pa- 
rameter the contribution (|29p can be comparable to lower 
orders in the expansion in l/k F L and b) the order pa- 
rameter A(e) has a maximum at the Fermi energy (e = 0) 
and decreases on an energy scale e ~ Ao as one moves 
away from the Fermi level. One can also show that meso- 
scopic corrections given by Eqs. (|26|) . (|27[) and (|29|) al- 
ways enhance A(0) as compared to the bulk value Ao. A 
couple remarks are in order: a) the energy dependence 
of the gap is universal in the sense that it does not de- 
pend on specific grain details, b) the matrix elements 
I(e, e') play a crucial role, e.g. they are responsible for 
most of the deviation from the bulk limit. Finally we 
briefly address the interplay of mesoscopic fluctuations 
and parity effects (see [28| for a more detailed account). 
The Matveev-Larkin (ML) parity parameter A p [23j . a 
experimentally accessible observable, accounts for even- 
odd asymmetries in ultrasmall superconductors. While 
the ML parameter coincides with the standard supercon- 
ducting gap in the bulk limit, in [23j j it was found that 
its leading finite size correction is given by 



A 



p = j^2N+1 — -^{E2N 



'->2N+2 I 



A(0)--, (30) 



where En is the ground state energy for a superconduct- 
ing grain with N electrons. 

We see that these corrections to the BCS mean-field 
approximation are comparable to mesoscopic fluctuations 
but have an opposite sign. For Al it seems that meso- 
scopic corrections are larger than those coming from (1301) . 



C. 2d chaotic cavities 

In this section we study a 2d superconducting chaotic 
grain of area A, perimeter C, and linear size L = \f~A. 
Our starting point is the gap equation ((3|) together 
with the semiclassical expressions for the spectral den- 
sity, Eqs. © and ©, and the matrix elements, 7(e,e'), 
Eqs. ([T0HI4|) . namely 



A 






1 


~ A 


C 




2nk F 


.4 


+IL 





L 

k F A 



C 



Si(4k F L) 



Ci 



. f4(e-e')k F L 



Ci 



V2(e-e') 



(F 



CF 



(31) 



where C « 0.339... and Si(x), Ci(x) are the sine and 
cosine integral functions, respectively. For I ^> L, the 
chaotic classical dynamics leads to a universal form for 
the function IL(u>), 



4 sin(k F luj) 
n 'H = TT ■ 



(32) 



As in the 3d case, the energy dependence of matrix el- 
ements implies that the equations to be solved for the 



gap are of integral type, and that the gap itself is en- 
ergy dependent. However, unlike the 3d case, we have 
logarithmic corrections coming from the contribution of 
the matrix elements. Based on the expansion in pow- 
ers of l/k F L of the spectral density and J(e, e') [see also 
Eqs. (|A.34I) and (IA.44[) ] we propose for a 2d chaotic grain 
the expansion 



A(e) = A [i + Z^+ZW+tt- 1 /^) 



(33) 



Following the same steps to solve the gap equation as in 
the 3d case, we get to leading order, 



(log) 



C log 2k F L 
2nk F A 



(34) 



Similar logarithmic corrections to residual interactions in 
2d chaotic quantum dots in the Coulomb Blockade regime 
were reported in Ref. 

The next order correction is given by 



a/ (1) = (c'±r 



c 



2k F A 



9(0), 



(35) 



with (— ) for Dirichlet and (+) for Neumann boundary 
conditions, respectively. The truncated spectral density 
g(0) is defined as in the 3d case, with semiclassical am- 
plitudes corresponding to 2d systems. 

Finally, the energy dependent correction to the gap in 
2d chaotic grains, f^ 2 \e) is given by the same function 
(I2TJ1) as in 3d grains. 

We note that a) in 2d the leading finite size contri- 
bution comes from the interaction matrix elements, not 
from the spectral density, b) finite size effects are stronger 
than in 3d and the leading correction does not vanish for 
any boundary condition, c) since effectively there are two 
expansion parameters l/k F L < 1- assuring the validity 
of the semiclassical approximation- and S/Aq < 1 - in 
order to apply the BCS formalism- it can happen that 
in a certain range of parameters the contribution /^(e) 
is dominant. 

In Fig. [3]we plot the gap as a function of the energy in 
units of the bulk gap Ao for Al grains (k F 17.5nm _1 , 
A w 0.18, and 8 7279/A meV where N is the number 
of particles), of different sizes L. Note the single peak 
at the Fermi energy. For the smallest grains the leading 
contribution is (e). This is yet another indication that 
the matrix elements play a dominant role in the finite size 
effects in superconducting metallic grains. 



IV. ENHANCEMENT OF 
SUPERCONDUCTIVITY IN NANOGRAINS: 
IDEAL VERSUS REAL GRAINS 

According to the findings of previous sections the su- 
perconducting gap is an oscillating function of the sys- 
tem size and the number of electrons inside the grain. 
Even for grains with N ~ 10 4 — 10 5 electrons consider- 
able deviations from the bulk limit are observed. For a 
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1.3 




FIG. 3. Superconducting order parameter A(e), Eq. (|33|l . 
in units of the bulk gap Ao for 2d chaotic Al grains (hp = 
17.5nm _1 ,5 = 7279/AT, A « 0.24meV) as a function of the 
energy e with respect to the Fermi level e = 0. Different 
curves correspond to grain sizes (top to bottom) and bound- 
ary conditions: L = 6nm, k F L = 105,<5/Ao = 0.77) (Dirich- 
let and Neumann boundary conditions), L = 8nm, = 
140,5/Ao = 0.32 (Dirichlet), and L = 10 nm, k F L = 
175, 8/ Ao — 0.08 (Dirichlet). The leading contribution comes 
from the energy dependent matrix elements I(e, e). 



fixed grain size, the deviations from the bulk limit are the 
larger the more symmetric the grain is. This is a typi- 
cal shell effect similar to that found in other fermionic 
systems, such as nuclei and atomic clusters [jj . These 
shell effects have their origin in the geometrical symme- 
tries of the grain. Symmetries induce degeneracies in the 
spectrum and, consequently, stronger fluctuations in the 
spectral density. The superconducting gap is enhanced 
if the Fermi energy is in a region of level bunching (large 
spectral density). Likewise, if the Fermi energy is close 
to a shell closure (small spectral density) the supercon- 
ducting gap will be much smaller than in the bulk limit. 

Therefore, thanks to shell effects, one can adjust the 
gap value by adding or removing few electrons in such a 
way that the Fermi energy moves into a region of high 
or low spectral density. In fact, shell effects in metal- 
lic grains of different geometries have recently attracted 
considerable attention [3 [lj, Q3, GJ, H Wj§ • A super- 
conducting spherical shell and a rectangular grain were 
studied numerically in Ref. [Tol A similar analysis was 
carried out in Ref. [lj for a nanowire. A qualitative analy- 
sis of a spherical superconductor was reported in Ref. Il7l . 

Discrepancies with experiments are expected because 
factors such as decoherence, deformations of the shape 
of the grain, and surface vibrational modes are not taken 
into account in the theoretical analysis. In this section we 
discuss the impact of small deformations of the grain and 
of decoherence effects that shorten the coherence length. 
We will see that weakly deformed grains can be modeled 
as symmetric ones but with an effective coherence length 



that incorporates the details of the deformation. The 
semiclassical formalism utilized in this paper is especially 
suited to tackle this problem. 

A. Superconductivity and shell effects 

We study the dependence of the gap on the number 
of electrons N inside the grain and compare the gap be- 
tween two grains with slightly different degree of symme- 
try. We focus on 3d rectangular grains where deviations 
from the bulk results are expected to be larger. In this 
case the chemical potential can be computed exactly as 
a function of N by inverting the relation 

\N = j\{e)de (36) 

where u{e) is the spectral density. 
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3.1xl0 4 3.2xl0 4 n 3.3xl0 4 3.4xl0 4 

FIG. 4. The superconducting gap A in units of Ao ~ 
2.286meV, as a function of the particle number N for a cu- 
bic (circles), of side L, and a parallelepiped-shaped (1.0288 : 
0.8909 : 1.0911) (squares) grain. Fluctuations are on average 
stronger in the cubic grain due to its larger symmetry. The 
parameters utilized are A = 0.3, eo = 32meV, e F ~ 11.85eV, 
k F = 26nm~ 1 . The energy gap was obtained by solving 
Eq. (|16|) with the semiclassical expression of the spectral den- 
sity given by by Eqs. (|A.7IA.8IA.10I) and a single particle 
coherence length I ~ 12L. 

As it is shown in Fig. 1 matrix elements does not 
affect the gap oscillations. Therefore we can solve the 
gap equation Q following the steps of section IIIII with 
the spectral density given by Eqs. (|A.7IA.8IA.10I) and 
/ ~ 1/V. The spectral density depends on the cutoff, 
namely, on the number of periodic orbits taken into ac- 
count. This cutoff is set by the single-particle coherence 
length I. Here we take I ~ 12 L where L is the length 
of the longest side of the parallelepiped and study the 
differences between a cubic and a rectangular grain. The 
cutoff is chosen to be much larger than the system size 
in order to observe fluctuations but considerably smaller 
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than the superconducting coherence length £ in order to 
accommodate other effects (see below) that might reduce 
the typical single-particle coherence length in realistic 
nanograms. We study a range of N such that the BCS 
theory is still applicable but deviations from the bulk 
limit are still important. 

In Fig. 2] we plot A, from Eq. (J3j) , as a function of N 
for a cube an a parallelepiped with aspect ratio 1.028 : 
0.89 : 1.091. For both settings we observe strong fluctua- 
tions with respect to the bulk value. The fluctuations are 
clearly stronger in the cubic case since the grain symme- 
try is larger. We also observe that a slight modification of 
the grain size (or equivalently N) can result in substan- 
tial changes of the gap. The observed differences between 
the cube and the parallelepiped are due to the different 
symmetry of these grains. In the cube the overall sym- 
metry factor in the spectral density is oc N 1 / 2 . The par- 
allelepiped has only two symmetry axis and therefore the 
symmetry factor ~ TV 1 / 3 . 

In addition to the fluctuations due to periodic orbits, 
we also expect smooth corrections to the bulk limit due 
to the surface and perimeter term of the spectral density. 
These corrections will be clearly observed as the coher- 
ence length is shortened and the contribution of periodic 
orbits is therefore suppressed. 
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FIG. 5. Superconducting gap A for a cubic grain (volume 
iV/181nm 3 ) for different single particle coherence lengths 
I = 2.25L, I = 6L, I = WL in units of A » 0.228meV 
as a function of the number of particles N. The param- 
eters utilized are A = 0.3, en = 32meV, £f ~ 5.05eV, 
k F = 18nm _1 . The energy gap was obtained by solving 
Eq. (|16p with the semiclassical expression of the spectral den- 
sity given by Eq.©. As the coherence length is reduced less 
periodic orbits contribute to the spectral density and fluctu- 
ations are smaller. Fluctuations are strongly suppressed for 
coherence lengths I < 2L. In this limit the gap is still smaller 
than Ao as a consequence of the surface and curvature terms 
in Eq.Cn>]) . 



B. Finite size effects in real small grains 

Highly symmetric shapes are hard to produce in the 
laboratory. It is thus natural to investigate to what 
extent small deformations from a perfect cubic shape 
weaken the finite size effects described in previous sec- 
tions. For applications it is also important to understand 
the dependence of the results on the single particle co- 
herence length I. In order to study this dependence, we 
assume that the superconducting coherence length £ is 
the largest length scale in the system. This is the most 
interesting region because in the opposite case I » £ the 
results for the gap (|2~T1) are to a great extent independent 
of I, By contrast, in the limit £ 3> I, the cutoff (j2TJ)) in- 
duced by £ has little effect as the contribution of periodic 
orbits L p > £ is already strongly suppressed by the cutoff 
induced by I. If I ~ £ both cutoffs must be taken into 
account. 

We now address these two related issues. We note that 
not only the effect of a finite coherence length I but also 
small deviations from symmetric shapes can be included 
in our analytical expressions for the gap by adding an 
additional cutoff V (besides Eq. (|20|) ) which suppresses 
the contribution of periodic orbits longer than T>. The 
details of V depend strongly on the source of decoherence 
or the type of weak deformation. Indeed, in certain cases 
T> may modify not only the amplitude but also the phase 
of the contribution of the periodic orbit to the trace for- 
mula used to compute the spectral density. For instance 
the effect of small multipolar corrections to an otherwise 
spherical grain [38] is modelled by adding an additional 
D cutoff in term of a Fresnel integral that smoothly mod- 
ulates the amplitude and phase of the periodic orbits of 
the ideal spherical grain. 

If the deformation is in the form of small, non over- 
lapping bumps, (39[, the cutoff is exponential and only 
affects the amplitude. The numerical value of the cutoff 
depends on the original grain and is directly related to 
the typical size of the bump. If the source of decoherence 
is due to finite temperature effects, [io}, T> - - — 



sinh(L p /Z) 

with I inversely proportional to the temperature. 

In Fig. [5] we show the effect of a finite coherence 
length I in the superconducting cubic grain investigated 
previously. The gap equation Eq. ([3]) was solved ex- 
actly with the semiclassical spectral density given by Eqs. 
(|A.7IA.12IA.10p and I = 1/V. For simplicity we use 
V = sin h( P /'/;) as a cutoff with I now the single parti- 
cle coherence length. This is enough for a qualitative 
description of the suppression of shell effects as a conse- 
quence of decoherence or geometrical deformations. 

The cutoff Eq. ([20|) . related to the superconducting co- 
herence length, does not affect the calculations as it is 
much longer (~ 1600 nm) than the ones employed in Fig. 
[5] Similar results are obtained if the analytical result 
dHJ) is utilized. 

As expected, the amplitude is reduced and the fine 
structure of the fluctuations is washed out as the coher- 
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ence length is shortened. We did not observe any gap 
oscillations with N for I > 2.5L. This can be regarded 
as an effective threshold for a future experimental verifi- 
cation of shell effects in superconductivity. Smooth non- 
oscillatory corrections depending on the S (or perimeter 
C in 2d) term in the spectral density are not affected by 
the coherence length and should be clearly observed in 
experiments. Note that A in Fig. [5] is, on average, below 
Ao even for the maximum N investigated. This is a di- 
rect consequence of the negative sign of the surface term 
in the spectral density for Dirichlet boundary conditions 
used in the numerical calculations (/^ in Eq. (ITT)) ). 



— is the mo- 



for 1/kpL <C 1, where kp — k(ep) 
mentum at the Fermi energy tp and L is the linear sys- 
tem size. The resulting semiclassical expansion will be 
organized in powers (possibly fractional) of the small pa- 
rameter 1/kpL. 

In order to observe deviations from the bulk limit, the 
single-particle coherence length must be larger than the 
system size, I > L. The time scale, r ss l/vp, associ- 
ated with I has a meaning of the lifetime of states near 
the Fermi energy. The condition I > L means that the 
Cooper pairs are composed of quasiparticles with a life- 
time longer than the flight time through the system. 



CONCLUSIONS 



1. Density of states 



We have determined the low energy excitation spec- 
trum, E — A(e) 2 + e 2 of small superconducting grains 
as a function of their size and shape by combining the 
BCS mean-field approach and semiclassical techniques. 
For chaotic grains the non-trivial mesoscopic corrections 
to the interaction matrix elements make them energy de- 
pendent, which, in turn, leads to a universal smooth en- 
ergy dependence (|29|) of the order parameter A(e), see 
Fig [3] In the integrable (symmetric) case we found that 
small changes in the number of electrons can substan- 
tially modify the superconducting gap, see e.g. Fig [4] 
Due to its potential relevance for experiments, we have 
investigated how these shell effects decrease (Fig. [5]) when 
the grain symmetry and/or the single-particle coherence 
length are reduced. 
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Appendix: Semiclassical approximation for the 
density of states and the interaction matrix 
elements. 



Semiclassical techniques such us periodic orbit theory 
[l| are not a common tool in the study of supercon- 
ductivity however they are a key ingredient in our an- 
alytical treatment. In order to solve the gap equation 
Eq. (|3]) we first need a closed expression for the spec- 
tral density and the interaction matrix elements I(e, e'). 
In this Appendix we describe in detail how these quan- 
tities are computed using a semiclassical approximation 



We start with the analysis of the density of states. 
The semiclassical expression for v(e) for a given grain 
geometry is already known in the literature [lj, 



~ ^tf(O) [l + g(e)+W)} 



(A.l) 



The spectral density gets both monotonous g(e) and os- 
cillating g(e) corrections. The monotonous correction at 
the Fermi energy is given by the Weyl expansion. 



9(0) 



Stt i 2C 
-ik F V " r Wv 



■ 2k F A 



3d 
2d 



(A.2) 



for Dirichlet (— ) or Neumann (+) boundary conditions. 
In Eq. (|A.2j) . S is the surface area of the 3d cavity, C is 
its mean curvature, while C is the perimeter in the 2d 
case. 

The oscillatory contribution to the density of states 
is sensitive to the nature of the classical motion. For 
a system whose classical counterpart is fully chaotic it is 
give n to the leading order by the Gutzwiller trace formula 



2tt 



-kpL p 



3d 



k F A L-i 



p A p e 



«(0 = ft ■ 

i[k P L p +^] e i^k P L p 2d, 

where we used k(e') ~ kp + e'kp/2ep. The summation is 
over a set of classical periodic orbits (p) of lengths L p < I. 
Only orbits shorter than the quantum coherence length 
I of the single-particle problem are included. The am- 
plitude A p increases with the degree of symmetry of the 
cavity [l[ (see below). In the chaotic case A p = A p (ep) 
is given by 



A p (e F ) = 



|det (M p - 1 



,1/2 > 



(A.4) 



with the monodromy matrix M p taking into account the 
linearized classical dynamics around the periodic orbit. 
The classical flow also determines [l| the topological in- 
dex P p in Eq. 
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Note that Eqs (|A.3[) and (|A.4[) indicate that the scaling 
of g in terms of the small parameter 



C = l/k F L 



m( e ') « 



C 2 3d, 
C 2d. 



(A.5) 



(A.6) 



Finally, for we have periodic orbits labeled by a single 
integer n 



47Taj 



■ Y!l^ cos ( k FL l n + ^k F L\^j 3d 
f^En^s (k F Li+^kpLi) 2d 

(A.12) 

with lengths L l n = 2na,i. The dependence on £ in this 



Rectangular grain 



3d 
2d. 



(A.13) 



Consider a rectangular box of sides a, with i = 1, . . . , d 
in d dimensions. For these systems the sum over periodic 
orbits is exact and given by [2, E| , 



S(0 



f5 (3) (e') 

-iEiE^sScO 3d, 



(A.7) 



2d. 



Here g^ is a sum over families of periodic orbits. Each 
family is parametrized by three (non simultaneously zero) 
integers n = (m, n?, n.3) 



i 



ff (3) (e') = E MkpLa + - kpLji) 



(A.8) 



where = 2-^/ a^n^ + a\n 2 + a^rig is the length of an 
orbit in the family and jo(x) = smx/x is the spherical 
Bessel function. We see that 



3 (3) oc C- 



(A.9) 



(2) 

In the same spirit, g\ ■ is written as a sum over families 
of periodic orbits parallel to the plane defined by sides 
cii , Oj . In this case the families are labeled by two integers 
n = (m, n 2 ) and 



(A.10) 



^ EUo Jb + 2d 
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k F LV + ^k F L l i> 

p n 2ep n 



3d 



where LV = 2 



a i n i + a 'j n 2 i s the length the orbit 
(711,712) and Jo is a Bessel function. Using the asymp- 
totic expression for J , we find that this contribution 
scales with £ as, 



C 3/2 3d 



€V)°c Ul/2 2d . 



(A.11) 



It is important to note that depending on the classical 
dynamics and the spatial dimensionality there are differ- 
ent types of scaling with £. The amplitude of the spectral 
fluctuations increases with the degree of symmetry of the 
cavity. It is maximal in spherical cavities and minimal in 
cavities with no symmetry axis [![. The latter typically 
includes chaotic cavities, namely, cavities such that the 
motion of the classical counterpart is chaotic. 

This relation between symmetry and fluctuations can 
be understood as follows. In grains with one or several 
symmetry axis there exist periodic orbits of the same 
length. As a result of taking all these degenerate orbits 
into account, the amplitude of the spectral density is en- 
hanced by a factor C^ 1 ! 2 for each symmetry axis |4lll42|. 
For instance, a spherical cavity has three symmetry axis 
so the symmetry factor is proportional to Q~ z / 2 1. 
Periodic orbits in chaotic cavities are not in general de- 
generate and the symmetry factor is therefore equal to 
one. For the range of sizes L ~ 5 — lOnm studied in this 
paper the difference between a chaotic and an integrable 
grain can be orders of magnitude. 



2. Interaction matrix elements 

a. Semiclassical approximation to the average density 

Unlike the case of the density of states, there is no 
general semiclassical theory for quantities, such as the 
interaction matrix element J(e, e'), involving the spatial 
integration of more than two eigenfunctions in clean sys- 
tems. For integrable systems the ergodic condition, 



A 




(A.14) 



with 51 = V or A in 3d and 2d respectively, is typically 
not met due to the existence of constants of motion. The 
constraints imposed by conservation laws effectively lo- 
calize the eigenfunctions in a smaller region of the avail- 
able phase space. 

On the other hand, for chaotic systems Eq. (|A.14|) is 
well justified as a result of the quantum ergodicity theo- 
rem [43| . The vast majority of the eigenfunctions spread 
almost uniformly over the whole volume (area) due to 
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the lack of constants of motion besides the energy. If the 
position r is far enough from the boundaries, we have 



l^(r)| 2 = ^(1 + 0(0) 



(A.15) 



for almost all states close to the Fermi energy. In or- 
der to evaluate explicitly deviations from Eq. (|A.14[) . we 
propose the replacement 



(A.K 



The average is over a small window of states around 
e n . The width of this window is controlled by an energy 
scale H/t related to the single-particle coherence length 
I ps vpT. This averaging procedure is justified since eigen- 
functions of classicall y ch aotic systems have well defined 
statistical properties |45| . 

The above average is exactly given by 



9(e) 
1 

7i"5(e) 



ffl(e')3G(f, r, e' 



iO + )de' 



where G{r,r',z) is the Green function of the non- 
interacting system at complex energy z, w(x) is a nor- 
malized window function of width H/t centered around 
x = 0, and g(e) is the density of states smoothed by w(x). 
Next, we express the Green function as, 



G = G° + G. 
where G° is given by the free propagator 



G°{r,f',e+iQ + ) = 



7) I 



2-nh 2 \r—f'\ 



3d 

r>\) 2d, 



(A.18) 



(A.19) 



and Hq is the Hankel function. The corresponding con- 
tribution to the average intensity, obtained by taking the 
limit r — > f 1 of the imaginary part of G° on lA.19( is then 
spatially uniform and given by 



(l^1| 2 )° = i. 



(A.20) 



The effect of such so-called zero-length paths joining f 
with r in zero time is then to produce a constant back- 
ground independent of the position (see, for example 
34]). This result should not come as a surprise, as zero- 
length paths are responsible for the leading order terms 
in the Weyl expansion of the density of states. 

In the semiclassical approach [29| the other part of 
the Green function, G, is expressed in terms of non-zero 
paths 7 going from r to r in a finite time r 7 as, 



G(f,r,e)=^£> 7 e 



i( kp- L-y-\- 2 ^ kF L-y~\-0-y 



(A.21) 



This contribution is responsible of the typical spatial os- 
cillations of the average intensity. The classical prop- 
erties of each trajectory are encoded in its topological 




phase f3 7 (equal to ir/4 times the number of conjugate 
points reached by the trajectory) and the smooth func- 
tion D 



(A.22) 



Here qi and q'j are local coordinates transverse to the tra- 
jectory 7 at points r, respectively, and f , and £ 7 (r, f ) is 
its length. In 3d we have two perpendicular components, 
while in 2d there is only one. 

After substitution of Eq. (|A.21|) into Eq. (|A.17|) . the 
integration over energies can be carried out explicitly pro- 
vided that, in consistency with the stationary phase ap- 
proximation used to derive the semiclassical Green func- 
tion, all smooth functions of the energy are evaluated 
at ep. The resulting Fourier transform of the window 
function acts as a cut-off for the sum. We finally obtain, 
after using the expression for the density of states and 
factorizing the Thomas-Fermi density, 



(Mr)\ 2 ) e = 



2 . _ 1 l + R(r,e) 



(A.23) 



In both 3d and 2d, R(r, e) is simply obtained form the 
Green function as a sum over classical paths 7(7^ = 7 
starting and ending at point r with finite lengths L 1 {r) = 
i 7 < I and actions S 1 {r) — hk(e)L^ [47j 



R(r, e) = ^D 1 cos ( k F L 1 + -^—kpL-y + /3 7 ) . (A. 24) 



Inspection of Eq. (|A.22I) shows that R scales as 

3d 
2d ' 

Furthermore, the normalization condition implies 




(A.25) 



i J R(r,e)df=g(e)+g(e). (A.26) 



Eq. (|A.26[) can also be used as the definition of the density 
of states without the Thomas-Fermi contribution. 

The separation between smooth, g(e) ~ g(0), and os- 
cillatory terms g(e) in Eq. (IA.26|) is as follows. Smooth 
contributions come from trajectories starting and ending 
at r after hitting the boundary only once L 7 < L. On 
the other hand, trajectories hitting the boundary more 
than once will have in general L p > L, and their contri- 
bution to the spatial integral can be evaluated using the 
stationary phase approximation to give g(e). 

Using Eqs. (|A.2IA.6IA.25|) and (|A~26l) the interaction 
matrix elements have the following semiclassical expan- 



sion, 



-2 2 
n 1 



3d, 



*[l+J( eF ,e,e')- T s • 
A [l + I{e F ,e,t l )~- 1 ^ + ~gi{e F )] 2d 



(A.27) 
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where 



R(r,e)R(f,e')dr. 



(A.28) 



b. Evaluation of I(e, e) 

As we will see, I(e, e') is a smooth function of both 
e and e'; it does not oscillate as rapidly as e lS / h where 
S is the classical action. The key point in carrying out 
the spatial integration in Eq. (|A.28j) is the separation 
of R = R shmt + R lons into short and long classical tra- 
jectories. A similar separation leads to the smooth and 
oscillatory contributions to the density of states discussed 
in the previous section. In other words, our approach to 
evaluating J(e, e') is similar to the Weyl expansion for the 
density of states. 

To calculate i? short , we note that in the regime I > L, 
the short trajectories of length L 1 < L are insensitive 
to the smoothing, and hence their contribution to the 
imaginary part of the Green function in (|A.17[) can be 
pulled out from the energy integration. This means that 
i? short is simply proportional to the imaginary part of the 
Green function associated with the short paths. 

Following Balian and Bloch [4l[ the basic idea of the 
subsequent calculation is that the boundary of the grain 
can be locally approximated as a plane in 3d and a 
straight line in 2d provided that the observation point 
is close enough to it. Within this approximation, the 
exact Green function representing a single reflection off 
an infinite wall can be calculated using the method of 
images. For the 3d case Eq. (|A.22[) is quantum mechan- 
ically exact and gives the same result. 

Following this idea, we construct the Green function 
for an infinite straight boundary by means of the method 
of images to obtain 



G 



ihort 



(r, 7, e + i0+) = ±G°(r, T(f),e + i0 + ) (A.29) 



where the action of the linear operator T is to map the 
position r* into its image point on the other side of the 
boundary. The plus and minus sign give Neumann and 
Dirichlet boundary conditions respectively. 

It is easy to see that when r — > r 1 , the function G short 
depends only on the distance between r and the bound- 
ary, which we denote by x. Since this distance is still of 
the order of the system linear size, it is possible to per- 
form the energy average in Eq. ([A. 171) with the Green 
function given by Eq. (|A~29|) . As a result, 



R 



short (z* 



e)=± 



{sin 2k(e)x 
2k(e)x 
Jo(2fc(e 



):/;) 



3d, 
2d, 



(A.30) 



After R shmt is inserted into Eq. (|A~28)) . the integral 
along directions parallel to the plane simply yields a fac- 
tor of S in 3d and C in 2d. The integration in the per- 
pendicular direction is naturally truncated at the system 



linear size L. In 3d, using f Q = — J^°, we obtain 



S TrSMin [k(e),k(e')} 



8k 2 F LV 



AV k(e)k(e') 



which, as expected, is a smooth function of e and e'. The 
second term in this expression was previously obtained 
in Ref. [2l| via a slightly different method which misses 
the first term of the right hand side of Eq. (IA.31I 

A similar analysis in 2d is more subtle due to diver- 
gence of the integration in the direction perpendicular to 
the boundary. However, there is a natural upper limit for 
this integration given by the linear system size L. Upon 
using k(e) ~ kp(l + e/2ep) and introducing the scaled 
perpendicular distance to the boundary y — 2k F x, 



rshort 
l 1d 



(e,e')= 



C 



2k F A 



2k F L 



Jo 



1 + £> 



Jo 



1+ 2^' y 



(A.31) 
dy. 



Employing the asymptotic expression for the Bessel func- 
tions, we find 



C 



c + - 



2k F A 

2kpL sin2y + cos2(e- e')y/e F 



(A.32) 



dy 



valid for kpL ^> 1 
^^(^-0.850 



"short 
2d 



case, /; 

e — e'). This implies 
perconducting gap is 
order in £. 

The integrals in Eq. 
the sine-integral (Si) 
Our final result is 



In Eq. (|A.32j) the constant C = 
. We see that, contrary to the 3d 
depends on energy (through the difference 
that in 2d chaotic systems the su- 
energy dependent even to leading 



(|A.32[) can be expressed in terms of 
and cosine-integral (Ci) functions. 



/ rshort /V ^ 



ttS 
ik F V 



\ rshort/, , ,A C 

hd ( £ F,e- e J - 7^4 



a 



Si(4fc F L) 



+ 



2ixk F A 



Ci ( 4(6 "^ )fcFL ) - Ci( 



3d, 



2d 



(A.33) 



with C = C- Si(2)/7r = 0.339 . . .. Thus for fixed e and 
e', 7 short scales with C = l/k F L < 1 as follows 



3d, 

l!^ rt ( eF ,e-e')c<C + &'Cic-gC 2d, 



7|f rt (e F )ocC + 6C 2 

rsho 
l 2d 



(A.34) 



where b and b' are constants independent of the system 
size. Note the non-algebraic dependence on ( in the 2d 
case. The constant b turns out to be much smaller than 
all other second order contributions to the gap, and will 
be dropped from now on. 

Now we focus on the contribution of long paths, i? long , 
to the spatial integral (IA.28|) . We use the expression 
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for R as a sum over classical closed paths 7(F) starting 
and ending at r with length L 1 (r) . Now we impose the 
condition 



L 7 (r) » L, 



(A.35) 



expressing the fact that the paths are long, namely, they 
hit the boundary several times. As is standard in these 
cases, we evaluate the smooth functions D 1 in i? long at 
the Fermi energy and expand fe(e) ~ kp + k(e) 2 /2kp to 
get 

7,7' 

(A.36) 

The phases involved in the spatial integration are (we do 
not include topological indexes for simplicity) 



$± 7 ,( e , e ',r) = fc F (L 7 ±L 7 -) 

+ |^(L 7 e±L Y e'), 



(A.37) 



where L 7 = L 7 (r) is the length of the trajectory 7. 

In chaotic systems different trajectories, in general, 
will have lengths differing by at least L (see, however 
48]). This means that since the first term in Eq. (|A.37j) 
scales as 1/C 3> 1, the integral over r in $Z"„/(e, e', r) and 
$ 77 (e, e',r) can be evaluated by the stationary phase 
method. Within this approximation, oscillatory integrals 
of the form 



f{x)e iXh ^dx 



are given to leading order in 1/ A by 



f(x)e lXh ^dx~f{x*) 



iXh(x* ) 



y/2TT\h"(x*y 



where h'(x*) — 0, and therefore each spatial integration 
in Eq. (|A.36|) yields an extra factor cx l/C 1 / 2 . Com- 
bining this with the prefactors (IA.25|1 . we find that the 
contribution of pairs 7^7' (the so-called non-diagonal 



contribution) is of order 



fiong/ n I C 5 ^ 2 3d, 
WfeOocj^ 2d . 



(A.38) 



On the other hand, terms that involve <£>" 7 (e, e',r) do 
not oscillate rapidly, because in this case the highly os- 
cillatory terms in the phase cancel each other leaving the 
second term in Eq. (IA.37[) which scales as £ and not as 

i/C 



$ 7)7 (e,e',r) = ^ ' (e-e). 



(A.39) 



This contribution involves coherent double sums over 
classical trajectories and is usually referred to as the di- 
agonal contribution, /j° ng . Taking 7 = 7' in Eq. (|A.36[) 



we easily find 

rlong / / \ 

7 d R fa (e^ ,e- e ) 



y^ J D 2 cos$^ 7 (e,e',r)dr, (A.40) 



which can be cast in a very compact form by introducing 
the purely classical function 

Ui(w) = J cos wk F L^(r)dr, (A.41) 



as follows 

Flong 



pong/ A _ I V V € F 



3d 
2d. 



(A.42) 



Keeping also in mind the ^-dependence of the coeffi- 
cients Z? 7 , we have 



/-long 



(e F ,e- e) oc 



C 2 3d, 
C 2d. 



(A.43) 



Equations (|A.33IA.38|) and (|A.42|) complete the eval- 
uation of I. Restricting ourselves to the first two orders 
in £ (C and Cl°gC m t ne 2d case), we finally obtain 

/(e,e')= (A.44) 



n 2 S 2 , /long/', , J\ 



3d, 



2d. 



Equations (IA.44I) together with the definitions (|A.33[) 
and (|A.42[) allow for the calculation of interaction ma- 
trix elements in 3d and 2d chaotic grains. In general, the 
explicit evaluation of II; (w) requires the precise knowl- 
edge of all classical paths up to lengths of the order of 
the single-particle coherence length / that have a cross- 
ing at r for every point inside the cavity. However, if I is 
large enough compared to L (in practice I ~ 5L suffices) 
ergodic arguments can be invoked and a closed expres- 
sion for the interaction matrix elements can be found. In 
situations when I ~ L one must carry out the explicit 
system-dependent calculation. 

Classical ergodicity of chaotic systems can be formu- 
lated in various ways [4{J, and we are going to give only 
a brief sketch of its consequences here. The main mecha- 
nism behind universality in the quantum mechanical de- 
scription of classically chaotic systems, resides in the be- 
havior of typical (in the sense of measure theory) classi- 
cal trajectories. By definition, a typical trajectory of a 
chaotic system will explore in an uniform way the avail- 
able phase space, thus implying the equivalence between 
temporal and microcanonical averages. 

This uniformity extends, in a non-trivial way, to the 
periodic orbits as well. The key concept here is the clas- 
sical probability of return, defined as 



P(x ,t,e) 



1 



Z(e,t) 



S(x -£(x ,t))S(H(x )-e) (A.45) 
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where xq = (rb,po) is a point in phase space mapped at 
time t into x(xq,£) by the solution of the classical equa- 
tions of motion. Clearly, the function 8(xq — x(xo,t)) 
is non-zero only when the classical flow maps an initial 
point into itself after a time t and plays the role of a 
probability of classical return. Moreover in case we want 
to select a fixed energy we use an extra condition given 
by the value of the Hamiltonian function along the tra- 
jectory. Finally, the probability must be normalized such 
that 

J P(x ,t,e)dx = 1 (A.46) 

thus fixing Z(t,e). The key observation here is that, by 
definition, the set of points where P(xo,t,e) is different 
from zero, belongs to periodic orbits with period t. Al- 
though the original ergodicity criteria were given in terms 
of typical trajectories, the theory of dynamical systems 
provides a strictly equivalent definition of ergodicity in 
terms of periodic orbits, 

P(xq, t, e) — > const. for t -> oo, (A. 47) 

That means that not only typical trajectories, but also 
periodic orbits uniformly fill the available phase space. 
We remark that the left hand side of this equation, a 
set of delta peaks at the periods of the classical periodic 
orbits, must be understood in the sense of distributions, 



namely, both sides are assumed to be integrated over a 
smooth function of time and phase-space position. 

In order to make contact with the coordinate represen- 
tation used so far, we use the uniformity of periodic orbits 
in phase space expressed by Eq. (|A.47I) , and integrate out 
the momentum. This integral can be exactly calculated 
[50j . It involves a Jacobian of the form dp(t)d(fo), which 
is indeed proportional to the semiclassical prefactors -D 7 . 
In summary, in the present context classical ergodicity 
leads to the following sum rule 50] for classical closed 
orbits, 

J2D? l 6(l-L, ( (r)) = \Y Z (A ' 48) 

7 I k F A M ' 

As was mentioned previously integration over lengths up 
to I on both sides with a smooth weight function is also 
assumed. Using this result and noting that the right hand 
side of Eq. (|A.48|) is independent of the position r, we get 

{4rr 2 sin wk F l 
y si nX* 2d ; ( a - 49 ) 

In the ergodic regime, I 3> L, these results enable us to 
evaluate explicitly the energy dependence of the interac- 
tion matrix elements in chaotic cavities. 
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